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Abstract 

We design an arbitrary high-order accurate nodal discontinuous Galerkin spectral element approximation for the 
nonlinear two dimensional shallow water equations with non-constant, possibly discontinuous, bathymetry on un¬ 
structured, possibly curved, quadrilateral meshes. The scheme is derived from an equivalent flux differencing 
formulation of the split form of the equations. We prove that this discretisation exactly preserves the local mass 
and momentum. Furthermore, combined with a special numerical interface flux function, the method exactly pre¬ 
serves the mathematical entropy, which is the total energy for the shallow water equations. By adding a specific 
form of interface dissipation to the baseline entropy conserving scheme we create a provably entropy stable scheme. 
That is, the numerical scheme discretely satisfies the second law of thermodynamics. Finally, with a particular dis¬ 
cretisation of the bathymetry source term we prove that the numerical approximation is well-balanced. We provide 
numerical examples that verify the theoretical findings and furthermore provide an application of the scheme for a 
partial break of a curved dam test problem. 

Keywords: split form shallow water equations, discontinuous Galerkin spectral element method, 
summation-by-parts, entropy stability, well-balanced, discontinuous bathymetry 


1. Introduction 


Fluid flows in lakes, rivers, and near coastlines are of interest in oceanography and climate modeling. For such 
flows the vertical scales of motion are much smaller than the horizontal scales. From this and the assumption of 
hydrostatic balance [T], the Euler equations can be simplified to the shallow water equations. If the fluid flows over 
a non-constant bottom topography the shallow water equations may be written as a hyperbolic system of balance 
laws 

wt + fx + gy = s. ( 1 . 1 ) 


It is well-known that solutions of the balance laws (1.11 may develop discontinuities in finite time, independent of 


the smoothness of the initial data. Hence, we consider solutions of the balance laws (1.11 in a weak sense that are 


well-defined provided the source term s remains uniformly bounded, i.e., weak solutions of (1.11 are well-defined 
under the assumption that the function used to model the bottom topography is in the space 1F''^(M), see e.g. [5]. 


The design of numerical methods to approximate o is driven by the need for stable, accurate and robust 
behaviors. For instance, the preservation of steady-state solutions is critical in problems with non-constant bottom 
topographies. Preserving steady solutions discretely is particularly troublesome for discontinuous bottom topogra¬ 
phies where special discretisations of the source term are required, e.g. mm- One steady-state constraint for the 
shallow water equations is the “lake at rest” condition II Eli], since the relevant waves in a flow can be viewed as 
small perturbations of the lake at rest, see [S]. A good numerical method for the shallow water equations should 
accurately capture both steady states and their small perturbations (quasi-steady flows) so as to diminish the ap¬ 
pearance of unphysical waves with magnitude proportional to the mesh size (a so-called “numerical storm” 0 ), that 
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are normally present for numerical schemes that cannot preserve the “lake at rest” condition. A numerical method 
that exactly preserves the “lake at rest” steady state property is said to be well-balanced, see e.g. [SJISlIilH]. 

Another critical requirement of the numerics is robustness and the ability of the method to remain stable 
and accurate, particularly the removal of aliasing errors that can drive nonlinear instabilities, and maintenance of 
stability even if discontinuities develop. Recent work has appeared on the use of high-order discontinuous Galerkin 
(DG) approximations to create robust numerical methods for the solution of systems of conservation laws, e.g., 
[El 13 HD]- These robust high-order DG methods may be derived from the perspective of mathematical entropy 
conservation, e.g. isnmiig, or reformulations of the PDE into a split formulation to maintain conservation, e.g. 
[3 EH]- The motivations behind the two approaches are similar [Id] . 

The split form of an equation is usually found by averaging its conservative form and non-conservative advective 
form. This is problematic as it is not obvious that discretisations of the split form remain conservative, yet 
conservation is desired for the numerical solution to exhibit correct shock speeds. Recent success has been had 
using diagonal norm summation-by-parts (SBP) finite difference operators to discretise the spatial derivatives in 
the split formulation of the equations [3 m EH Eg. Fisher et al. E3] show that split form operators derived 
from SBP derivative matrices are consistent and conservative in the Lax-Wendroff sense. There is now a known link 
between SBP finite difference operators and the discontinuous Galerkin spectral element approximation with Gauss- 
Lobatto points, e.g. m- This link was used in [3 to derive an entropy conserving discontinuous Galerkin spectral 
element method (DGSEM) for the one dimensional shallow water equations. This paper exploits the links further 
and extends in a non-trivial way the previous work found in [5] to multiple dimensions and possible discontinuous 
bottom topographies. 

In this paper we present an entropy stable, high-order discontinuous Galerkin spectral element approximation 
for the shallow water equations with a discontinuous bottom topography for unstructured and curved quadrilateral 
grids. The DGSEM is naturally discontinuous at element boundaries, so we ensure high-order (spectral) accuracy 
by placing element boundaries at discontinuities in the bottom topography. The ability to do so allows one to 
model realistic bottom topographies appearing in oceanography. The scheme presented here is also well-balanced, 
an attribute difficult to guarantee in curvilinear coordinates. We find that the numerical satisfaction of the metric 
identities El] (referred to in m as the geometric conservation law) is critical to prove that the baseline scheme 
remains entropy conservative and well-balanced on arbitrary meshes. 

Our approach is to use results of Fisher m and Fisher and Garpenter ESI to derive an entropy conserving 
approximation, and from that an entropy stable one, which is possible because we can reformulate the spectral 
element approximation of the split form of the shallow water equations into an equivalent flux differencing structure. 
We use the flux differencing reformulation to prove the underlying properties of the entropy stable DGSEM, as well 
as to highlight how an existing DGSEM code can be altered to incorporate entropy stability. 

The paper is organised as follows: in Sec. we begin with a brief description of the entropy analysis of the 
two dimensional shallow water equations. We outline the discontinuous Galerkin spectral element method with 
the summation-by-parts (SBP) property in Sec. This section also introduces the important reformulation of 
the DGSEM into an equivalent flux differencing framework, which is the critical equivalence that allows us to use 
existing theory. We provide in Sec. |4.1| a discretisation of the two dimensional shallow water equations using the 
flux differencing formulation that is conservative and entropy conservative on curvilinear meshes. We also provide a 
detailed proof that the approximation remains well-balanced. Then in Sec. |4. 2 [ additional dissipation is added to the 
scheme to ensure that the approximation remains valid for flow regimes that may contain shocks. Numerical results 
in Sec. demonstrate and underline our theoretical findings. Our conclusions are presented in Sec. Finally, 
[Appendix Gj provides algorithms and implementation details of how a standard DGSEM code can be altered to 
incorporate the newly proposed entropy stable fluxes. 


2. Shallow water equations 

We begin with the balance law form of the two-dimensional shallow water equations 

ht {hu)x + {hv)y = 0, 

{hu)t + {hu^ +g /i^/ 2 ) 3 , -k {huv)y = -ghb^, (2.1) 

{hv)t + {huv)x + {hv"^ + gh?/2)y = -ghby, 
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which includes the continuity equation and the momentum balances. The quantity h = h(x,y,t) denotes the water 
height measured from the bottom topography b — b{x,y) with the total height given hy H = h + b. Additionally 
the constant g is the gravitational acceleration. The fluid velocities are given by u = u{x,y,t) in the a;—direction 
and V = v{x,y,t) in the j/—direction. The SW model (2.11 is compactly written as a system of balance laws 0 
with w = {h, hu, hv)'^, the fluxes 


/ = {hu, hu^ + g y?/2, huv)'^, g = {hv, huv, hv^ + g 


( 2 . 2 ) 


and the source term s = (0, —ghb^, —ghby)"^ 


Since the system (2.1) is nonlinear, we must define in what sense our numerical approximation will be stable. 
The extension of the usual stability for linear problems is the so-called entropy stability m, where a (generalized, 
mathematical) entropy function rather than the norm of the solution is non-increasing in time. To this end we 
impose the entropy condition as an additional admissibility criterion on the system. 

The entropy condition states that for smooth solutions the entropy of the system is conserved and for discon¬ 
tinuous solutions the entropy decays. Numerical approximations that satisfy the entropy condition discretely are 


referred to as entropy stable. The balance law (1.1) for the shallow water equations is equipped with a convex 


mathematical entropy function e = e{w) in the form of the total energy [5] 

e := (^2 + ^, 2 ^ 


(2.3) 


To develop the conservation law for the entropy function e(w) we define the set of entropy variables, q = 
(li, 92 , 93 )^, by 

de 


<?2 := ^— = u , 
OW2 

de 

<?3 := ^— = V. 
0 W 3 


(2.4) 


To determine the entropy fluxes, iF{w) and G(w), we use two compatibility relations that must hold between them, 
the entropy variables and the physical fluxes m 


^w — q fwt Gw — q Qw- 


(2.5) 


We pre-multiply the balance law (|l.l|) with the entropy variables (|2.4) and apply the conditions (2.5) to And the 

( 2 . 6 ) 


conservation law for the entropy function 

Ct + J-x Gy = 

Explicitly, the entropy fluxes of the shallow water equations are 


T := + huv^) + g{hu{h + b)), 

G ■= + ^w^z)) -I- g{hv{h + b)). 


In the presence of discontinuities the entropy conservation law (2.6) becomes the entropy inequality 

et + iFx Gy ^ 0. 


(2.7) 


( 2 . 8 ) 


We then define an entropy conserving approximation of the nonlinear shallow water equations, (1.1), to be one that 


discretely satisfies (2.6) and entropy stable if (2.8) is satisfied 


Finally, we note that systems of balance laws have important steady state configurations where the flux and 
source terms are in balance. For the shallow water equations, one such steady state solution is the “lake at rest” 
condition defined by 

h + b = const, 

n (2.9) 

u = V = U. 

A numerical method that preserves the “lake at rest” state is said to be well-balanced. If a method is not well- 
balanced spurious waves on the order of the mesh size truncation error can be generated and pollute the approxi¬ 
mation. 
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3. Nodal discontinuous Galerkin spectral element method 

The entropy stable method that we propose is a form of nodal Discontinuous Galerkin Spectral Element Method 
(DGSEM). In this section, we introduce the basic construction of the DGSEM on curvilinear quadrilateral grids. 
A more complete discussion can be found in m- Also in this section, we provide details about the relationship of 
DG methods to summation-by-parts (SBP) operators, which allows us to write the approximation as an equivalent 
sub-cell flux differencing formulation (PDF) [TU]. The PDF is useful for theoretical purposes, particularly the proof 
of local conservation and satisfaction of the Lax-Wendroff condition. Therefore, the DGSEM has local and global 
conservation. Also, special choices of the flux functions in the PDF generate a DG discretisation of a split form of 
the original PDE. Additionally, if the underlying flux functions in the PDF are two-point entropy conserving fluxes 
then the PDF remains high-order and entropy conservative m- Because of the equivalence between the PDF and 
the DGSEM with the SBP property this automatically generates an entropy conservative DGSEM. 


3.1. Conservation law in curvilinear coordinates 


The DGSEM approximates the system of conservation laws (1.1) defined on a domain 17 on an unstructured 
mesh of quadrilateral elements. To simplify the discussion, we work with the components of the system (2.1) written 
as 

(.Wk)t +ifkiw))^ +{gkiw))y = Skiw), fc= 1,2,3, (3.1) 

where wi = h, W 2 = hu and W 3 = hv. 


We decompose 12 into non-overlapping quadrilateral elements G and for computational efficiency map each 
element to the computational reference element E = [—1,1]^. A commonly used transformation between the 
reference square and an arbitrary curve-sided quadrilateral element is transflnite interpolation with linear blending 
HZ]. The mapping between the coordinates of the reference square (^, rj) and the physical coordinates x = {x, y) is 


v) = \ [(1 - e)r4(r7) + (1 + 0 H 3 ) + (1 - ^)f i ( e ) + (1 + v)r3{0] 

- J [(1 - ?){(! - ^)ri(-l) + (1 + r?)f3(-l)} (3.2) 

+ (l + a{(l-ry)fi(l) + (l + ry)f3(l)}]. 


where we assume that each element is bounded by four curves P^, j = 1, 2, 3,4. 

Under the transformation the conservation law on 12 remains a conservation law on E. We transform the two 
dimensional balance law (3.1) from physical space to the reference space by rewriting derivatives using the chain 
rules 

dw dw dw dr] 

dx d^ dx dr] dx’ , , 

dw dw d^ dw dt] 

dy d^ dy dy dy' 

for some w{x,y). Applying the chain rule (3.3) to gives us the transformed conservation law in reference space 

Jiwk)t + ifk)i + igk)'n = 0, (3.4) 


where the element mapping (3.2) is used to compute the Jacobian, and the contravariant flux components fk, gk, 
k = 1,2,3 according to 

J = x^yr,- Xyy^, 

fk{w) = yy fk{w) - Xygk{w), (3.5) 

gk{w) = fk{w) + x^ gk{w), 

for fc = 1, 2, 3. 
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3.2. Polynomial approximation on the reference element 

The DGSEM is constructed by approximating the variables Wk and contravariant fluxes fk, Pk in the reference 
space by polynomials of degree N in each direction. We use a nodal form of the interpolation with nodes defined at 
the Legendre-Gauss-Lobatto (LGL) points {Ci}^o {?7j}^o i"^ the reference square E = [—1,1]^. The Lagrange 
basis functions for the interpolant are 




N 


n 

2 =0, 




J = 0, 


.N. 


(3.6) 


and satisfy the cardinal property 




(3.7) 


where Sij denotes Kronecker’s symbol with Sij = 1 for i = j and Sij = 0 for i + f. We write the element-wise 
polynomial approximation (e.g for the components of w) as 


Wk{x,y,t)\^ = Wk{x{^,r]),y{f,r]),t) 


N N 

Wkit V, t) := E E (A m A iv), k = 1 , 2 , 3, 
2=0 3—0 


(3.8) 


where j=o the time dependent nodal degrees of freedom. The nodally represented polynomial 

( |3.8[ ) is equivalent to an orthogonal Legendre polynomial expansion used in a modal spectral method, but is more 
convenient to use in the approximation of nonlinear equations like the shallow water equations. 


We use the idea of collocation throughout this work to approximate quantities derived from the Wk. For instance 
the velocity u is approximated by a polynomial of degree N in each direction (3.8 1 as well, where its nodal values 
are computed as 

z,j = 0,...,W (3.9) 




Wf 




This collocation strategy also applies to the contravariant fluxes, where we interpolate the metric terms at the same 
nodes. For instance 


K’' = vAA, Vj) FkiW^A - xAA, Vi) GkiW^A, i,j = 0,...,N, A = 1,2,3 


(3.10) 


where the three flux components for the shallow water equations are defined in (2.2 1 . Similarly, the Jacobian of the 
transformation is approximated by the polynomial of degree N with nodal values 

For spectral approximations, the derivative is approximated elementwise directly from the derivative of the 
polynomial approximation, e.g.. 


r) N N „ 

^ i=0 j=0 ^ 

where k = 1,2,3. We introduce the polynomial derivative matrix D with entries 

dij 


D — 

- df, 


i, j = 0,..., IV, 


(3.11) 


(3.12) 




which is used to calculate the derivative with respect to f. Since we use the same polynomial ansatz in ^ and rj 
direction, the derivative operator is identical in each direction. 


We can reuse the ID operator (3.121 for the 2D scheme if we store the individual slices in a 2D array. Entries 


in one column refer to the nodal values at constant 77 , entries in one row share the same ^ 


\WkA 


= -l,r?=-l 


{WkA 


=-1.j7=+1N 


'W 


0.0 


W, 


0,N ' 


w. := 


(3.13) 




=+l,r)=-l 


{WkA 


=+l,r)=+l 




Af.O 


w, 


N,N 


By storing the nodal values of variables, fluxes and metric terms in such 2D arrays, we can multiply the ID derivative 
operator D from the left to represent taking the f derivative at each constant 77. If we multiply with from the 
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right, we obtain the 77 derivative. This notation allows us to write the scheme in a compact matrix-vector notation 
while using unmodified ID operators. 


With notation (3.13), the nodal values of the derivative of a nodal approximation w are given by 


or in index notation 


where j = 0,..., iV. 


(Wfe)^ = DWfe and (Wfe),, = WfcD^, 


N 


N 


= E iWl’% = DuWi’^ 


1=0 


1=0 


(3.14) 

(3.15) 


We demand that the discretisation preserves free-stream solutions, i.e. constant solutions of the balance law 


(3.4) should remain constant for all times. A necessary and sufficient condition for constant state preservation is 

(3.16) 


that the metric identities 


d d ^^2 s' 

a -I- a = (1, 

0 ^ or] 


are satisfied, where the volume weighted contravariant basis vectors, = 1,2 are 


, Ja^ = . 


(3.17) 


The metric identities are not automatically satisfied for the discretisation, which, using the array notation (3.131 
for the nodal values of the metric terms and (3.14) for the derivative can be expressed as 


Dyr, - = 0 

—Dx,, -I- = 0. 


(3.18) 


Kopriva |20j proved that free-stream preservation is guaranteed for the DGSEM when the linear blending formula 


(3.2 1 is used and boundaries of the quadrilateral elements are approximated by polynomials with an order equal to 


(or lower than) the polynomial order of the approximate solution. Thus, we use an isoparametric approximation 
in which each boundary curve rj(^), G [—1,1] of an element G is approximated by a polynomial of order N. We 


use the same Lagrange basis functions (3.6) to approximate the boundary curves 


N 




(3.19) 


3=0 


where due to their robust interpolation properties m, the nodes {Cil^o typically chosen to be the Chebyshev- 


Gauss-Lobatto or Legendre-Gauss-Lobatto nodes. The polynomial boundary curve approximations (3.19) are used 


to construct the mapping (3.2) for each element. As the mapping is a polynomial in ^ and rj, the derivatives 


necessary to obtain the metric terms and the normal vectors are computed directly in the discrete derivative sense 


(3.14). Further details of the isoparametric polynomial approximation of boundary curves can be found in nil HO]. 


3.3. Discontinuous Galerkin spectral element method (DGSEM) 

Following a standard approach we derive a nodal discontinuous Galerkin scheme, e.g. [T7| or [01]. We will 
discuss the specific form of the approximation of the source term in Sec. ^ Omitting the source term, the nodal 
discontinuous Galerkin method in weak form of the transformed conservation law (3.4) reads 


f J {Wk)t^d(,dri - f Fkip^d^drj- f Gkifn d(,dr] = - (f (f^(W^,W ),Gl(W'^,W )]-hipdS, 

JE,N Je,N Je,N JdE,N ^ ' 

(3.20) 

for k = 1, 2, 3 and in the equivalent strong form m 


[ {J{Wk)t + {Fk)i + {GkY) ^d^d7] = -(f (F:{W+,W-)-Fk,Gl{W+,W-)-Gk) -nyrds, (3.21) 

Je,N JdE,N ^ ' 
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for A: = 1,2, 3, where is a polynomial of degree N or less in each space direction. We also use the notation ^ 
to represent Legendre-Gauss-Lobatto Quadrature, which is the tensor product of one space dimension quadrature 


1 


-l,iV 


N 


j =0 


(3.22) 


with being the LGL quadrature weights. The quadrature is exact if the integrand / is a polynomial of 

degree 2N — 1 or less. The numerical fluxes and G^, in the normal direction, couple neighbouring elements. We 
indicate this coupling by the dependence on the “outer” and “inner” values W^,W~ along the normal vector n. 


The boundary integrals in (3.201 and (3.211 describe the integration along the four edges of the element E. 


The two forms (3.201 and (3.211 of the DGSEM are algebraically equivalent because the quadrature satisfies a 


summation-by-parts rule m- In one space dimension, and for any two polynomials U (^) and V (^), exactness or 
the quadrature implies that 


In other words. 


JV .1 .1 N 

^ = / UV^d^ = UVt, - / U^Vd^ = UVt, - ^ U',V,u;,. 

j =0 “'-1 “'-1 j =0 

[ UV^di = t/F|y - [ U^VdC 
J-1,N J-1,N 


(3.23) 


(3.24) 


The result extends to two and three space dimensions [^ . 

The formal statements of the DGSEM, ( |3.20[ ) and (3.21) can be reduced to a pointwise form m, which in turn 
can be represented in a matrix form where the nodal values are represented as arrays. To get equations for the 
nodal degrees of freedom, we take (/? = Then, for example. 


lE,N 


J{Wk)ty:^d^dr] = 


N 

E 

n,m—0 


UJnU:rnJ"'^{W^nd^i^n) iVrn) = ^,j = 0,1,2, 


,iV, 


(3.25) 


for fc = 1, 2, 3. 


as 


We can represent each component (3.25) in terms of the matrix-array notation introduced in (3.13) and (3.14) 

(3.26) 


MJo(Wfc),M, 

where M is the diagonal matrix of the quadrature weights (the mass matrix), 

M :=diag(wo,---,WAr), (3.27) 

and we introduce the notation of a component-wise Hadamard product of two matrices with the same dimension 

A o B = C, with dj = Qij bij, i,j = 0,..., N. (3.28) 

Similarly, the integral of the k — th component of the ^ contravariant flux in strong form is 


[ {Fk),ipd^dr] = W f E (?") iVm) = DiA , 

Jen V;=o / Vi=o / 


so we see that 


lE,N 


{Fk)^(pd^dri -> MDFfcM. 


On the other hand, for the weak formulation, 


[ Fkip^d^dr]= ^ U}nU}mFj!’'^i'i{^n)ij{Vm) = 

n,m—0 \n—0 / 


(3.29) 


(3.30) 


(3.31) 
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so 


f E,N 


Fw^d^drj —> D^MFfcM, 


(3.32) 


with a similar result for the ry direction flux. 


The boundary quadratures for the weak form, (3.20), are 


/ {Fk,Gl) ■ hd^(fdT] = f F^ ipdr] - [ ipdr] 

JdE,N J-1,N J-1,N 

+ [" Gl{^,l)d^- r GU^,-l)^dt 

J-l.N J-1,N 

Each term reduces to pointwise values in the same way. The first term, for instance, is 

r*! rl 


(3.33) 


[ (1, ry) (/3dry = / F^ {l,r)) £,{1) (rj) dr] = uJj{F^) 

J-1,N J-1,N 


N,j 


(3.34) 


We represent the numerical fluxes and consistent to the 2D array notation (3.13) this time so that the only 
non-zero entries correspond to the respective interfaces, i.e. 


/(^*)0,0 


(^*)0.Arx 



0 

0 


((Gir° 0 • 

■ 0 (G*) 

K ■■= 

0 

0 

, G*:= 

0 • 

* 

. . . O 


0,N ' 


N,N 


(3.35) 


so that F^ only appears at ^ = ±1 and appears at 7y = ±1 on the reference element E. 
To write the approximations in a compact form, we define the matrix operators 

D := —M“^D^M (scaled derivative matrix), 

S := diag [—,0,...,0,- — ) (surface matrix), 

Vwo UJnJ 


(3.36) 


with the derivative matrix D defined in (3.12). Then we can rewrite the equations for the nodal degrees of freedom 

(3.37) 

J o (Wfe)t -k D Ffc + Gfe = S (FI - Ffe) + (G^ - G^) S, (3.38) 


of the DGSEM approximations (3.20) and (3.21) in the algebraically equivalent forms 

J o (Wfe)t-(-D Ffc-k Gfe = S n + G^ S, 

and 


for k = 1,2,3. We note that to obtain the results (3.37) and ( |3.38[ ) we multiplied by the inverse of M on the left 
and right. 

3.4- An equivalent subcell flux differencing form 

The most important property the DGSEM operators constructed with the LGL quadrature nodes have is that 


they are summation-by-parts (SBP) operators for all polynomial orders, (3.24). This property can be represented 
in the form of SBP-SAT finite difference operators [H], and we collect relevant results here in Lemma 

Lemma 1 (SBP-Properties). Let the matrix Q MD, which represents the mass weighed derivative, as seen in 
{3.30). The matrix Q has the SBP-property 


Q = B := dm5(-l,0,...,0,l) (SBP). 

Furthermore, the SBP-property can be used to obtain alternative expressions for the derivative matrix 


D = M"iQ = M-i(B - Q^) = -S - 
= (M-iQ)^ = -S - QM-b 


(3.39) 



















So, the derivative matrix D of the weak DG formulation 


satisfies the relations 


D = 

(3.40) 

D = -S-kD, 

= -S-kD^. 

(3.41) 


Proof See, for example, [H [13 EEl IMl HI] • ® 

Remark 1. Although the finite difference and spectral element approximations differ (e.g. the spectral element 
approximate solution is known everywhere, including in between the nodes), the fact that the nodal equations can 
be written in the same form will allow us to simply use results proved for SBP finite difference methods as needed. 

The LGL-based DGSEM operators listed in Lemma are in the sub-class of SBP operators with diagonal 
norm matrix M. For this class of diagonal norm SBP operators, Fisher and Garpenter m proved an astounding 
relationship: Such operators can always be re-written into an algebraically equivalent subcell finite volume type 
differencing formulation. As an example we rewrite the contravariant flux in the x—direction into a telescoping flux 
form (the contravariant flux in the ?/—direction has an analogous form): 


DFfc = M-iQFfe = M-iAFfc, 
where A is the A^ x TV -|- 1 differencing matrix 


(3.42) 


/-I 

0 

0 
0 
VO 


1 


0 


-1 1 


0 

0 

-1 

0 


0 0 \ 

0 0 

0 0 
1 0 
-1 ly 


(3.43) 


The new flux functions, denoted with an overbar, can be viewed as subcell ffrrxes on a complementary staggered 
subcell grid miin]. The contravariant flux functions on the complementary grid remain consistent and high-order 
when they are computed according to m Section 4.5 and Appendix A.3] 


pO,j _ poj 
~ -^k ^ 


N 2-1 


^k ~ 


m—i 

^k ~ k ^ 


^ ^ 2 TF™-^) {{ Ja} {{Jal]] , * = 1,..., TV, 


(3.44) 

for fc = 1,2,3, a fixed point j in the y—direction and for some two point, symmetric flux functions and G™^, 

0'g' ^ 

(3.45) 


The F^'^ and F^'^ are the typical contravariant fluxes. The metric terms are given in (3.171, and the arithmetic 
mean is defined as 

(3.46) 




For computational efficiency we generalise a previous result of Fisher and Garpenter da Fq. (3.13)] and rewrite 
the flux differencing volume term on curvilinear meshes. 

Proposition 1 (Flux Diffferencing with Metric Terms). One can use the structure of the SBP matrix Q and 


the fluxes on the complimentary grid to eliminate one of the sums in (3.441 to write the flux difference in the 
X— direction, M“^AFfc, in the indicial form 


pi+ij _ pi,3 

f k 

iOi 


1 

UJi 


N 

2Q,„ (F^ylF-^iF-^■) + GrV^F^MF™’^) {{ J , (3.47) 


m—0 
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for i = 0^..., N. The flux difference in the y-drection, A^M ^, can be expressed in a similar indicial form 


a 


ij' + l 


- 


, N 

= - E 2Q.™ (Fr\w^’\w^n , (3.48) 


m—O 


/or j = 0 ,..., N. 


Proof. The details to derive (3.471 and (3.481 are straightforward and can be found in Appendix A 


With the definition of the flux differencing formulation complete, we must select the specific form for the volume 
fluxes and Depending on the choice of the volume flux it is possible for the flux differencing scheme to 

recover discretisations of alternative split forms of the PDE. The split form of a PDE, often found by averaging 
the advective and conservative form of the equations, is known to have stabilisation properties for non-linear 
PDE discretisations [ISl Hi HI]. But, due to their form, it is often unclear if the approximation remains globally 
conservative in the sense of Lax-Wendroff. We see from the telescoping flux difference formulations (3.471 and 
(3.48) that it is trivial to show conservation of the scheme while maintaining the flexibility and positive stabilisation 
properties of the split form. 

In particular, a split form for the derivative of a product of two quantities is 


(a b)x = ^{ab)^ + ^{aa;b + a b ^), 


(3.49) 


To demonstrate the split form property of the flux differencing form we select the form for the first component of 
the volume flux, say to be 


-pvol 


1 




(3.50) 


We substitute (3.50) into the first term on the right of (|3.47 1 and after some manipulation obtain 


Af .. / N N 

* \m—0 m—0 


OJi 


m—0 


N 


(3.51) 


+ {JalY’^ Y. 


m—0 


which is the i—th. row for the discretisation of (3.49) 

1 


{JalWO, 


(D (Ja] o Wi) -k Ja] o DWi -k Wi o DJa() . 


(3.52) 


That is, the product of two averages in (3.47) creates a discretisation of the standard split form of a quadratic 
product. In a similar fashion the product of three averages creates a discretisation of the standard split form of a 
triple product proposed by Kennedy and Gruber m- This is a remarkable property of the flux differencing form 
(3.47). By inserting an arithmetic mean or products of arithmetic means into the flux differencing scheme creates 
a discrete version of a particular split form of the equation. Complete details and proofs of this property of (3.47) 
can be found in pS] . 

So, rewriting the volume contributions of the DGSEM into the flux differencing form grants us additional 
flexibility to construct an approximation that discretises alternative split forms of the PDE. By using the alternative 
split form of the shallow water equations it is possible to create an entropy conservative numerical approximation [5] . 
This gives us the motivation to select the internal volume fluxes in such a way that the total energy of the numerical 
scheme will be conserved discretely. We note that the only alteration needed to change a standard DGSEM code 
to an entropy stable one is to change the volume contributions to the flux differencing form and select appropriate 
volume fluxes, see [Appendix C| 


To summarize, we have rewritten the volume contributions of the standard strong form DGSEM (3.38) into an 
equivalent flux differencing framework 


J o (Wfc)i + M-i AFfc + Gfc A^M-i = S (F* - F^) + (G^ - G^) S, k = 1, 2,3, 


(3.53) 
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where the flux differencing is computed in each direction according to the form (3.47) or (3.481. The flux differencing 
representation guarantees that with a chosen pair of symmetric volume fluxes the approximation (3.531 remains 
high-order and conservative. Additionally, if the volume fluxes in (3.471 and the surface fluxes in (3.381 are carefully 
constructed then the approximation is also provably entropy conservative m- Note that one can independently 
select an entropy conservative volume flux that has a different form than an entropy conservative surface flux. 
We will show in the next section that the additional flexibility to select different fluxes in the volume and at the 
surface allow the construction of an entropy conservative approximation for the shallow water equations that is also 
well-balanced. 


4. Entropy stable DGSEM for the 2D shallow water equations 

In this section, we construct the volume and surface fluxes, following the ideas of Tadmor, e.g. min], that 
discretely conserve the entropy. We denote the entropy conservative approximation as the ECDGSEM. Entropy 
conservation is only valid, however, for smooth solutions and not discontinuous ones (shocks). In Sec. |4.2[ we 
use the entropy conserving scheme as the starting point and add dissipation in a controlled way to guarantee that 
entropy is always dissipated at shocks, resulting in an entropy stable approximation denoted by ESDGSEM. 


4-.1. EGDGSEM on curvilinear meshes 


The flux differencing formulation of the strong form DGSEM on curvilinear meshes (3.53) provides a compact 
notation for the ECDGSEM for which the previous results of Fisher and Carpenter [T3] for conservation and entropy 
conservation apply. To ensure that the approximation remains well-balanced we use the extra flexibility of the flux 
differencing form that allows us to select different entropy conservative volume and surface fluxes. The discretisation 
of the source term is also divided into volume and surface contributions. The surface parts depend on averages, 
{{•}}, and jumps, !•], across the interface. We denote the jumps for an arbitrary nodal quantity W consistent to 
notation (3.13): 


0,0 


[WIj := 


0 




0 




IWI, := 


'iwir 0 

0 


0 

0 


(4.1) 


where we distinguish between jumps in ^ direction, IWJ^, and ry direction, IW],,. The averages of a nodal quantity 
are defined analogously. Whereas the local average operators are symmetric and hence don’t prefer a direction, we 
define the local jumps according to the ^ and rj coordinate directions: at the ^ = — 1 and the rj = —1 interfaces we 
compute the local jumps as the “inner” value minus the “outer” value, whereas for the ^ = 1 and rj = 1 interfaces 
the local jumps are computed as the “outer” value minus the “inner” value. 

Theorem 1 (Gurvilinear EGDGSEM). The semi-discrete flux difference form of the two dimensional EGDGSEM 
for the shallow water equations on curvilinear grids 


J o (Wfe)t + M-i AFfe + GfcA^M-i = S (F*’“ - Ffe) + (G*’“ - Gfc) S H- source^, fc = 1, 2, 3, 


(4.2) 


where for the flux differencing components (|3.47|) and (|3.48|), we use the entropy conserving volume fluxes 




(4.3) 




{{hv}}^,(3,„,) {{w}} 
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in combination with the entropy conserving surface numerical fluxes 






V {W}{M}{M} 
{W}{M}{W} 


(4.4) 


where in {4-4) the {{■}} indicates the average of the two neighbouring states W~^ and W~ with the source term 

\ 


discretisation 

( 


0 


-f h o [y^ o Db + D(y„ o b) - yj o bD^ - (y^ o b)D^] + fM ^ (^y^ o {{h}}^ o [bj^j - f (^y^ o {{h}}^ o [b|„^ M ^ 
^-fh o [-x„ o Db - D(x^ o b) + xj o bD^ + (xj o b)D^] - (x,, o {{h}}^ o [bj^j + f (^xj o {{h}}^ o [bj^j 


has the following properties: 


(4.5) 


1.1 Discrete conservation of the mass and discrete conservation of the momentum if the bottom topography is 
constant. 

1.2 Discrete conservation of the total energy, which is an entropy function for the shallow water equations. Hence 
it preserves the entropy of the system. 

1.3 Discrete well-balanced property for arbitrary bottom topographies. 

Proof. We prove the result in three parts. 

Proof of Part o The discrete conservation of the numerical scheme follows directly from the telescoping flux 
differencing form of the approximation m- 

Proof of Part |1.2[ If the volume and surface fluxes are provably entropy conservative then the global entropy 
conservation of the approximation is retained by the flux differencing form m- The surface fluxes ( |4.4[ ) are known 
to be entropy conserving [2j. We demonstrate in Appendix B.l that the volume fluxes (4.31 are also entropy 
conservative. 

Proof of Part |1.3| The complete proof of well-balancedness for the curvilinear ECDGSEM is provided in [Appendix] 
|B.2| However, it is necessary to describe the construction of the bottom topography discretisations. For the two 
dimensional problem in general coordinates we require discrete approximations for b^ and by, or more compactly, 
V&. From (3.3) we know the explicit form of the gradient in computational coordinates is 

T 


JVh = 


db 


db 


db 


db 


^^dr] 




(4.6) 


(4.6) is 


We treat each piece of the source term as a quadratic split form (|3.49|). For example the first component from 

(4.7) 


bcr. - ^ 


1 f _ 1 f&a+„,»+igsv 


2 \ df 


df df J 2 \ dr] 


Vi 


dt] dr] 


We then approximate the bottom contribution from the discrete approximation of the quadratic split form (3.52). 
So, we have the approximations for the derivative of the bottom topography 


1 


)D(y^ o b) -b o Db -(- b o Dy,, - (y^ o b)D'^ - yj o bD^ -bo y^D^ 


(4.8) 


by ~ 2 o b) — o Db — bo Dx,, -|- (x^ o b)D^ -I- xj o bD^ -|- b o x^D^) . 
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The bottom topography discretisation could immediately be treated in a flux differencing way with the result (3.47). 
However, we know that the approximation satisfles the metric identities (3.16). So we cancel extraneous terms in 


(4.8) and obtain a more computationally efficient form of the source term volume contributions in (4.2) 


1 


- (D(y,, o b) + y,, o Db - (yj o b)D^ - y^ o bD"^) , 
by K, - (^—D(x,, o b) — X,, o Db + (x^ o b)D^ + xj o bD”^) 


(4.9) 


The surface contributions follow from a similar logic in that we require the derivative of b at the boundary (in this 
case the jump). 


Remark 2 (One-dimensional bottom formulation). In one space dimension we use the following approximation of 
the source term for an element G 


g hbx ^ gi^Db' 




^ Wo 






cat+i. 

UJn 


(4.10) 


Here is the average water height at the left interface node of element G and the right interface node of 

element G — 1, is the average water height at the right interface node of element G and the left interface 

node of element G -I- 1: 




(4.11) 


In the same fashion the jump in bottom topography is defined as 


mf^+i = - b%. 


(4.12) 


The quadrature weights are denoted by uiq and lon and the -|- 1 dimensional unit vectors are denoted by ei and 
e/v+i- 


4--2. Entropy stable DGSEM 


The ECDGSEM presented so far exactly conserves the discrete entropy. However, the solution of the shallow 
water equations may develop discontinuities (shocks) in finite time even for smooth initial data. We know in the 
presence of discontinuities that the conservation law for the entropy function (2.6) must be replaced by the entropy 
inequality (2.8) m- Thus, we must add numerical dissipation to the ECDGSEM so that the entropy function 
is guaranteed to be dissipated (or conserved for smooth well resolved solutions), thereby ensuring that a discrete 
version of the entropy inequality holds. 


A typical way to add dissipation in a discontinuous Galerkin approximation is via the definition of the numerical 
flux function. Most often, those numerical flux functions are (approximate) Riemann solvers that inherently create 
dissipation at shocks (or where the solution is otherwise underresolved). We follow this basic idea and add dissipation 
in the spirit of Riemann solvers at the element interfaces to add dissipation to the ECDGSEM. 


To derive the dissipative numerical flux, we note that the physical fluxes (2.2) have the associated flux Jacobians 


and 


= fw = 






(4.13) 


(4.14) 
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The eigenvalues of (4.131 are u + c, it, u — c and for (4.14) v + c, v, v — c with the wave speed c = '/gh- The 
matrices of eigenvectors of (4.13) and (4.14) are 


Rf = u + c 0 u — c 


(4.15) 


and 


respectively. 


Rg — 


1 

u 


1 

u 


(4.16) 


lU + c 0 V — cj 


The dissipation term will also require the entropy Jacobian matrix H = To obtain H, we now express the 
conservative variables w = (/i, hu, hv)'^ in terms of the entropy variables q = {g{h + b) — u, u)^ : 

Wi=-qi-b + ^{ql + ql), 

9 2g 

W 2 = -qiq 2 -bq 2 +—{ql +q 2 ql), (4.17) 

- ^93 + ^(^293 + qD- 


Differentiating (4.17) directly leads to the entropy Jacobian matrix 


H = - u gh + 


V 

uv 


uv gh + 


(4.18) 


With an appropriate scaling for the right eigenvectors we obtain the set of entropy scaled eigenvectors m 

H= (RT)(RT)'^, (4.19) 

which relates the right eigenvectors to the entropy Jacobian. For the scaling, we consider the matrix 

T = diag( v^, v^), (4.20) 

with scaling parameters on the diagonal only. We define Z = and have the identity in a new form 

H = RZR^. (4.21) 

For the eigenvectors of the / flux Jacobian, R^, we find 


Si = 


1 


S2 = h, S3 = 


1 


(4.22) 


A straightforward calculation shows that the same scaling works for the eigenvectors of the flux Jacobian in the 
y—direction as well. 

Now we have all the necessary components to define the entropy stable numerical flux functions. We subtract 
the dissipation term required for dissipation in the x—direction 


^ p*,ec_ |A^.| ZRj|g] 


1 , 


(4.23) 

(4.24) 


and the y—direction 

G*’^* = G*’“-iRg|Ag|ZR^|<fl, 

where Af and Ag are the diagonal matrices containing the eigenvalues previously computed. We use the arithmetic 
average values at an element interfaces to compute the right eigenvector, scaling, and eigenvalue matrices in (4.23) 
and (4.24). 
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It is important that the dissipation terms depend on the jumps of the entropy variables and not on the jump 
of the conserved quantities as would be common in standard Riemann solver-based numerical flux functions. If we 
compute the discrete entropy equation by multiplying the scheme with the entropy variables, we get contributions 
of the form — ^ |q] - R/ |A/| Z RJ |q] at each interface that are guaranteed to be negative due to the positivity of 
the matrix R/ | Ay | Z Ry . Writing in terms of the jumps in the entropy variables ensures that entropy is dissipated 
when the jump in entropy variables across interfaces is large (e.g. shocks) and is nearly preserved when the jumps 
are small for well resolved smooth solutions. 

We Anally present the main contribution of the present work, an an entropy stable DGSEM (ESDGSEM) for 
the shallow water equations. 

Theorem 2 (Gurvilinear ESDGSEM). The semi-discrete form of the two dimensional ESDGSEM formulation for 
the shallow water equations on curvilinear grids is given by 


Jo(Wfc)t + M-iAF, 


GfcA^M-i = S(F:’""-Ffc) 


fGr*-Gr 


S-l-sourcefc, A:= 1,2,3, 


(4.25) 


where the flux differencing components use the volume fluxes dot in combination with the entropy stable numerical 
fluxes (4.231 and (4.241 and the source term discretisation (4.5l. The approximation has the following properties: 


2.1 Discrete conservation of the mass and discrete conservation of the momentum if the bottom topography is 
constant. 

2.2 Discrete entropy stability. 

2.3 The well-balanced property for arbitrary bottom topographies. 


Proof. The ESDGSEM follows directly from the curvilinear EGDGSEM presented in Thm.j^ To guarantee entropy 
stability we replace the entropy conserving numerical fluxes (4.41 at element interfaces with the entropy stable 
numerical fluxes (4.231, (4.24). For the “lake at rest” initial conditions the jump in entropy variables is zero, |(f] = 0, 
so the additional dissipation term vanishes and does not affect the well-balanced property of the scheme. ■ 


5. Numerical results 


In this section, we use the entropy conserving (4.2) and entropy stable (H .25[) numerical schemes on several 


(|4:^ 

gTc 


test cases to numerically verify the theoretical findings from Thms. and To integrate the systems in time 
we use the five stage, fourth order Runge-Kutta time integrator of Garpenter and Kennedy m- First, to verify 
the convergence, conservation and well-balancedness of the approximations we use a structured curvilinear mesh 
depicted in Figure Elements are numbered by counting from left to right and bottom to top. We also simulate 
a bore-shear interaction as well as the numerical generation of potential vorticity generated from the passage of a 
non-uniform bore. After all theoretical findings are verified numerically, we present a simulated partial dam break 
from a parabolic dam with a discontinuous bottom topography in the downstream region of the flow. The partial 
dam break problem serves to exercise each component of the ESDGSEM approximation. 


5.1. Convergence 

We first test the convergence of the entropy conserving and entropy stable approximations with a smooth solution 
test problem. We use the method of manufactured solutions to create an analytic solution 


H(x, y, t) = h{x, y, t) -I- bi{x, y) = 8-\- cos{x) sin{y) cos(t), 
u{x,y,t) = 0.5, 
v{x,y,t) = 1 - 5 . 


The manufactured solution (5.1) introduces additional source terms to the equations of the form 


51 := Ht + u(H^ - (6i)x) + v{Hy - {bi)y), 

52 := uHt -\- - (^i)a:) + uv{Hy - {bi)y) -\- H^{H - bi), 

5 3 := vHt -I- uv{Hx - (bi'jx) + v^{Hy - {bi)y) -\- Hy{H - 6i), 


(5.1) 


(5.2) 
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y 0 



-0.5 


Figure 1: The curvilinear mesh used for verification of convergence, conservation, and well-balancedness. 


where derivatives regarding H and bi are computed analytically. We solve this problem on the domain [—1,1]^ with 
the smooth bottom topography 

bi{x, y) = 2 + O.Ssin {2Trx) + 0.5cos {2Try). (5.3) 

The gravitational constant is set to g = 1. 

We vary the polynomial degree on the mesh given in Fig. ^and observe exponential convergence up to = 16 
{N = 15 for ESDGSEM) for At = 1/2000 and N = 17 {N = 16 for ESDGSEM) for At = 1/4000, when the errors 
introduced by the time integrator become dominant. We present semi-log plots in Fig. [^for the entropy conserving 
scheme and Fig. [^for the entropy stable scheme. As previously observed, e.g. [EHHinilll], we find a suboptimal 
order of convergence for odd polynomial degree N for the purely entropy conserving scheme. However, both the 
EGDGSEM and ESDGSEM are spectrally accurate for smooth problems. 


5.2. Conservation of mass and momentum 


We first numerically verify that mass and momentum are conserved. Proper ty [T] l of Thm. fusing a constant 
bottom topography. Additionally, with the specific numerical volume fluxes (4.31 and surface flmces (4.4) the 
approximation will conserve the total energy (modulo dissipative effects of the time integrator). Also, for a non¬ 
constant bottom topography the momentum equations become balance laws, and we show that mass and entropy 
are still conserved discretely, even for discontinuous bottom topographies. 


The problem that we choose to model is a dam break on the domain H = [—1,1]^. The dam break is initialised 
along the vertical line a: = 0 on the curvilinear mesh shown in Fig. with periodic boundary conditions and a 
polynomial degree of A = 5. The gravitational constant is again set to g = 1. The dam break problem uses the 
initial conditions 


h{x,y,0) 


5 - b{x,y) 
4:-b{x,y) 


if a; < 0 
if a; > 0 ’ 


u{x,y,0) 


v{x,y,0) = 0. 


(5.4) 


5.2.1. Dam break over aflat bottom 

We demonstrate the entropy conservative properties of the EGDGSEM scheme. Property [Tl2 of Thm. We 
consider a flat bottom topography, 6=0, and solve the dam break on the curvilinear mesh Fig/^ The differences 
in mass, momentum, and total energy are listed in Table The error in the discrete energy reflects the dissipative 
influence of the time integrator. Otherwise, we see that the conservation in mass and momentum is on the order of 
machine precision for each time step value considered. By shrinking the time step we see that we can drive the error 
in the discrete total energy to machine precision, converging at the fourth order accuracy of the time integrator. 
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Figure 2: Semi-log plot shows the spectral convergence in space and fourth order accuracy in time for of the 
ECDGSEM scheme applied to a smooth solution. 



Figure 3: Semi-log plot shows the spectral convergence in space and fourth order accuracy in time for the ESDGSEM 
scheme applied to a smooth solution. 


5.2.2. Dam break over a diseontinuous bump 

Next we examine the conservative properties of the numerical scheme for a discontinuous bottom topography. 
Momentum will no longer be conserved. But the mass should be conserved and the error in the total energy should 
reduce as the time step is refined. These properties are demonstrated in the numerical test presented in Table 
For the discontinuous bottom topography we used 


b2{x,y) 


2 + 0.5 sin (27ra;) -I- 0.5 cos {2TTy ), if elementjjj = 6 
0, otherwise 


(5.5) 
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At 

A Mass 

A Moment umX 

AMomentumY 

AEnergy 

Temporal Order 

1/1000 

3.55E-14 

2.66E-13 

4.32E-17 

4.79E-08 

- 

1/2000 

2.49E-14 

2.66E-13 

9.95E-16 

3.01E-09 

3.99 

1/4000 

3.20E-14 

2.66E-13 

1.71E-15 

1.89E-10 

3.99 

1/8000 

3.20E-14 

2.66E-13 

1.46E-15 

1.18E-11 

4.00 


Table 1: Errors in the conserved quantities, mass, momentum, and total energy, at T = 1 over a constant bottom 
topography. In the energy conservation results we observe the temporal accuracy of the time integrator. 


which is the bottom topography (5.3) restricted to a single element. If we use the entropy conserving scheme without 
added stabilisation and periodic boundaries, we expect the entropy (total energy) to be conserved in the scheme. 
Table shows that the error in the total energy shrinks with the fourth order accuracy of the time integrator as 
the time step is refined. Also, we see that mass is conserved to machine precision for all time steps. 


At 

AMass 

AEnergy 

Temporal Order 

1/1000 

5.33E-14 

2.16E-08 

- 

1/2000 

1.78E-14 

1.35E-09 

4.00 

1/4000 

2.84E-14 

8.48E-11 

3.99 

1/8000 

3.55E-15 

5.32E-12 

3.99 


Table 2: Errors in the mass and total energy at T = 1 for the discontinuous bottom topography (5.51. In the energy 
conservation results we observe the temporal accuracy of the time integrator. 


5.3. Well-balancedness over a discontinuous bottom 


Next we demonstrate numerically that the curvilinear entropy conserving numerical scheme (4.2) is well- 
balanced, numerically demonstrating Property [^3 of Thm. We focus, particularly, on a discontinuous bottom 
topography. So, we configure a “lake at rest” test problem as in (2.9) 


h-\-b2{x,y) = 5, 
u = V = 0, 


(5.6) 


with the discontinuous bottom topography (5.5). The boundary conditions are set to be 
curvilinear mesh in Fig. Id on the domain = [—1,1]^ and vary the polynomial degree N. 
at At = 1/1000. Table 1^ shows that the L 2 -error of the approximate total water height, 
magnitude of round-off errors for both the ECDGSEM and the ESDGSEM. 


periodic. We use the 
The time step is fixed 
H = h b 2 , is of the 


N 

La-error EGDGSEM 

La-error ESDGSEM 

3 

8.84E-15 

5.37E-15 

4 

8.75E-15 

5.02E-15 

5 

1.85E-14 

1.55E-14 


Table 3: L 2 -error of the approximate total water height, H = h-\-b, to the “lake at rest” problem on the curvilinear 
mesh shown in Fig. |^at T = 1. 
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5-4-. Dam break over a discontinuous bump 

Next we compute the solution of a dam break problem using both the ECDGSEM and the ESDGSEM. This 
numerical example demonstrates that the entropy stable approximation removes spurious oscillations in the post¬ 
shock regions of the flow introduced by the entropy conservative approximation. For the numerical test we set 
inflow/outlow Dirichlet boundaries along the vertical lines x = 0 and x = 10 and periodic boundaries along the 
horizontal lines y = 0 and y = 10 on the domain fl = [0,10]^. We use a rectangular mesh with sizes varying from 
20 X 20, 40 X 40, 80 x 80 to 160 x 160 elements and polynomial degree of iV = 4. The gravitational constant is 
again set to y = 1. The initial conditions are 


/i(x,y,0) 


3.5 — &3(x, y) if X < 5 

2.5 — &3(x, y) if X > 5 ’ 


u{x,y,0) = v{x,y,0) = 0, 


(5.7) 


so we initialise the problem with a height discontinuity on the element interfaces at x = 5.0. As a bottom topography 
we use 


b3{x,y) 


2.0 — (x — 5)^ — {y — 5)^, if |x — 5| < 1 and |y — 5| < 1 
0, otherwise 


(5.8) 


which is a box with a smooth top that has its center at (5.0,5.0), side lengths of 2, and is initialised discontinuously 
along the edges of the box, which align with cell interfaces. 

The results shown in Fig. [^for the purely conserving scheme show that it produces severe ringing in the post¬ 
shock region of the approximation. The entropy stable approximation removes these spurious oscillations except 
near the discontinuity at the shock front. We show the computed entropy stable solution in Fig. where we present 
a grid refinement study for the entropy stable approximation. It is clear that the additional dissipation smoothes 
the spurious oscillations and the ESDGSEM provides a more physical solution to the dam break problem over the 
discontinuous bump. Nevertheless, we note that, although stable, the ESDGSEM is not completely oscillation-free 
0 . 



Figure 4: EGDGSEM, dam break over a discontinuous bump on 40 x 40 elements at T = 1 and CFL = 0.1. 


5.5. Non-linear breaking shallow waters waves 

This test problem simulates the interaction of a hydraulic bore in a shear flow. We also observe how the 
numerical scheme generates potential vorticity through the passage of a non-uniform bore iniisii. For smooth 
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(b) 40 X 40 


I 

I 


(a) 20 X 20 




I 

(c) 80 X 80 (d) 160 X 160 

Figure 5: Grid refinement study for the dam break problem modeled by the ESDGSEM at four grid resolutions 
with CFL = 0.1 and = 4. 






flows the potential vorticity 


A dv du 


(5.9) 


is a conserved quantity [SUES]. Hydraulic bores are discontinuities in the flow where energy is dissipated but 
mass and momentum are conserved. By design the dissipation in the ESDGSEM is generated proportional to the 
magnitude of the jump in the entropy variables which are large near discontinuities but spectrally small in smooth 
regions of the flow. Since the dissipation is only applied near bores, potential vorticity can be generated through 
non-uniform shallow water wave breaking. 


To see how the ESDGSEM generates potential vorticity we begin with an initial flow that has zero vorticity. 
We consider a flow where the bottom topography is zero, 6 = 0, with the initial linear gravity wave 


h{x, y) = 1 + Asin{ly) sin(fca;) 
u{x,y) = -^ sin(Zy) sin(fcx) 

(jj 

v{x,y) = —- cos{ly) cos{kx) 


(5.10) 


in the rectangular domain ft = [—0.5, 0.5]^. We set solid wall boundary conditions at y = ±0.5 and periodic 
boundary conditions in the a;—direction. The quantities in the initial conditions (5.10) are 


k = 27rm, I = (2n ± l)7r, = g{k^ + l'^), A = 0.1, 5=1, to = 2, n = 0. 


(5.11) 


We run three simulations for the initial conditions (5.10) each with 40,000 degrees of freedom: 


1. TV = 1 with a uniform 100 x 100 Gartesian mesh. 
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2. N = 3 with a uniform 50 x 50 Cartesian mesh, 

3. N = 7 with a uniform 25 x 25 Cartesian mesh. 

The first configuration makes the ESDGSEM a second order spatial approximation which can be directly compared 
to previous results, namely Fig. 21 in m and Fig. 4 in m- The other two configurations yield higher order spatial 
approximations to be used for comparison with the low-order computation. 

Due to the non-linearity in the shallow water equations, the higher amplitudes in the gravity waves begin to 
break at t « 0.5. The breaking waves create peaks at the crests and troughs of the free surface near to the walls, as 
shown in Fig. The breaking extends to the interior, and bores moving in the negative a;—direction are formed. 



(a) T = 0.0 


(b) T = 0.5 



Figure 6: Three dimensional visualization of the initial water height and the water height just before wave breaking 
occurs at t « 0.5. 

The bores are aligned in the y—direction with some curvature shown at T = 2 for iV = 1, = 3 and = 7 in the 

left part of Fig. We also present, in the right part of Fig. the point-wise potential vorticity at T = 2 computed 
using the local derivative to approximate the vorticity A. The computed water height and potential vorticity for 
A^ = 1 compares well with the results of Tassi et al. m- For the N = 3 and N = 7 computations we see there 
are some oscillations in the vicinity of the discontinuities and more resolution in the smooth parts of the flow. The 
generated potential vorticity follows the shape of the discontinuities in the solution. We again note that ESDGSEM 
adds dissipation to ensure entropy stability but is not guaranteed to be overshooot free. This explains the noise 
at the element boundaries in the potential vorticity plots particular for the N = 7 computation. Without any 
post-processing of the numerical solution visible numerical artifacts are aligned with the computational grid. 


5.6. Parabolic dam break 


In Secs. |5.1| - [5^ we have verified the theoretical properties of the EC and ESDGSEM. For the next example 
we combine each aspect of the numerical scheme and model the fluid flow from the partial break of a parabolic 
dam. First, we will demonstrate the well-balancedness of the ESDGSEM approximation on curvilinear grids, i.e. 
numerical verification of Property |^3 of Thm. To do so we consider the solution before the failure of the dam, 
which amounts to two “lake at rest” problems on the left and right sides of the dam. On the right side of the dam 
we also place a discontinuous bottom topography. Then, we allow the dam to fail and examine the flow. 


For both numerical tests in this section we use a domain O = [—5, 5]^, which is divided into 1600 quadrilateral 
elements. We model a parabolic dam placed near the center of the domain O with the curve 


1 


y 


1 

4' 


(5.12) 


Finally, for each configuration, we place a discontinuous bottom topography on the downstream side of the dam of 
the form 


b4{x,y) 


2.0-bln(a:-1.25) if a; > 2.25, 
0 if a; <2.25, 


(5.13) 
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5.6.1. ESDGSEM well-balancedness 


We first solve the “lake at rest” problem on either side of the dam before it fails. We will consider the water 
height to the left of the dam to be higher than the water on the right of the dam. We consider the initial conditions 


Hx,y,0) = 


1 

4’ 


10 - b^ix) ifx < 
5-64(2;) 


0 ) = i;(a;, y,0) = 0. 


(5.14) 


We set periodic boundary conditions for each lake individually. This test serves to demonstrate the well-balancedness 
of the ESDGSEM on a curvilinear mesh, including a discontinuous source term ( 5.13| ), Property]^ 3 of Thm|^ For 
the test problem we take the time step to be At = 1/5000 and integrate to a final time of T = 5. We present the 
L 2 error in the approximation of the constant total water height, H = h b, in Table We find that on either 
side of the parabolic dam the error in the computed water height is on the order of machine precision. 


N 

L 2 -error ESDGSEM (left) 

L 2 -error ESDGSEM (right) 

3 

5.52E-14 

4.21E-14 

4 

7.04E-14 

6.22E-13 

5 

1.82E-13 

1.05E-13 


Table 4: L 2 -error of the approximation of the total water height, H = h b, to the “lake at rest” problem for the 
parabolic dam at T = 5 for various values of N. The second column shows the error to the left of the curved dam 
where there is no bottom topography. The third column gives the error to the right of the curved dam where there 
is a discontinuous bottom given by (5.13). 


5.6.2. Partial dam break with discontinuous bottom topography 

The final demonstration considers the partial failure of a parabolic dam. The initial conditions are given by 


(5.14). The boundary conditions are periodic along the lines y = 5 and y = —5, Dirichlet along the lines x = 5 
and X = —5, and reflecting wall boundary states along the unbroken parts of the parabolic dam. We assume 
instantaneous failure of the portion of the dam in the region y G [—0.5, 0.5]. It is only in this region that the two 
states interact. 

First, we provide a visual grid convergence study for this complex test problem that has no analytical solution. 
In Fig. we provide the computed solution of the water height for three polynomial orders N = 3, N = 5, and 
N = 7, with At = 1/1500, integrated to a final time of T = 1.5. The overlay of quadrilaterals represents the spectral 
element mesh. We see from the numerical p-refinement study that the waves in the approximation are well-resolved 
for N = 5 and N = 7. 

From the grid convergence study we know that the computation is sufficiently well resolved with polynomial 
order A = 5 in each element. So with A = 5, we show the evolution of the water height of the partial dam 
break problem at times T = 0.0, T = 0.5, T = 1.0, and T = 1.5 in Fig. Again for this computation, we chose 
a time step of At = 1/1500. This numerical test combines each aspect of the ESDGSEM approximation, i.e., a 
discontinuous solution, curvilinear mesh and discontinuous bottom topography. The pseudocolor plots in Fig. 
show the propagation of eddies near the dam break. Lastly, we provide in Fig. [^a three dimensional visualization 
of the partial dam break where we can see on the downstream side of the dam the interaction of the flow with the 


discontinuous bottom topography (5.13). 


The numerical solution of the parabolic dam break problem demonstrates that the entropy stable numerical 
approximation can capture shock and rarefraction waves. However, the dissipation added to guarantee entropy 
stability is not designed to make the method overshoot free. Additional shock capturing techniques for DG type 
methods are necessary to remove the remaining oscillations, e.g. [50] . 

Lastly, we note that a standard DGSEM scheme is unstable when solving the partial dam break problem from 
a parabolic dam, even for very small time steps. 
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6. Conclusion 


In this work we developed a new high-order entropy conserving and entropy stable DGSEM discretisations for 
the two dimensional shallow water equations on general curvilinear meshes. To highlight the conservative property 
of the approximation on curvilinear meshes we reformulated the approximation into an equivalent flux differencing 
form. With this reformulation it is straightforward to demonstrate local conservation. Applying results of Fisher and 
Carpenter m, a careful choice of the numerical volume and surface fluxes leads to an entropy conservative scheme. 
The flux differencing form also provided additional flexibility to guarantee that the approximation remains well- 
balanced. For non-constant bottom topographies we found that the recovery of special steady-states of the shallow 
water equations depends on a special discretisation of the nonlinear source term. By considering a particular source 
term discretisation we maintained well-balancedness for discontinuous bottom topographies in general curvilinear 
coordinates. Finally, it is known that energy must be dissipated at shocks, but the entropy conserving scheme is 
dissipation free modulo any dissipative effects of the time integrator. The numerical solution can therefore capture 
shocks and rarefactions accurately but at the cost of significant post-shock oscillations. Thus, we provided entropy 
stable numerical fluxes to add dissipation to the scheme and control overshoots. Note that the dissipation added 
is merely the amount necessary to guarantee entropy stability and is not designed to make the method overshoot 
free. 

We provided six numerical examples to demonstrate and underline the theoretical findings. The simulation of the 
flow from a parabolic shaped partial dam break exercised each component of the novel FSDGSEM approximation. 


Appendix A. Proof of Prop. 


Proof. To show the flux difference formula (3.471 we first consider the single difference for the components i = 1 
and 1 = 2 and a fixed index for j. The argument presented readily extends to the other flux difference components 


in both the i and j directions. From the high-order flux extension on curvilinear grids (3.441 we have for i = 1 

N 


_ 

= E {{ J , (A.l) 


and for i = 2 


m—1 


N 1 


p2j _ 
- 




m=2 i=0 


We expand the second component (A.2) to find 

N 


p2j' _ 
— 


m ^2 

+ 2Qi^ {{ W^’^) {{Jai}} 




We subtract (A.l I from (A.31 and cancel like terms to determine 


_ pl.j _ _ 


-k 


N 


' (l,m) j' 

We know from the nearly skew-symmetric structure of the SBP matrix Q that 

~ Qoi = Qio, Qii = 0. 


j' 


(A.2) 


(A.3) 


(A.4) 


(A.5) 


We use the properties (A.51 and that the arithmetic mean and the volume flux, by assumption, are symmetric to 
collect the terms from (A.41 into a single sum 

_ _ AT 

= E {Fr\W^’^,W^’^) {{ Ja}{{ J . (A.6) 

m—0 
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We then generalise the calculation of the flux difference 


F, 


d+lj 


N 


771=0 


for i = 0,..., iV. We then premultiply by the inverse of M to obtain the desired flux differencing result (3.471 for 


M-iAF 


pi+l.i _ pi,j 
£fe 

UJi 


1 ^ 


771=0 


An identical strategy can be used in the j index direction to rewrite the flux difference in the ?/—direction, 
GfeA^M-\ in the similar indicial form (13.481). I 


Appendix B. Proof of Thm. 

Appendix B.l. Proof of Property^ 2 

Proof. In the flux differencing scheme we use different but consistent fluxes for the volume and interface contri¬ 


butions. The interface fluxes (4.4) are known to be entropy conservative [3]. We will demonstrate here that the 
volume fluxes are also entropy conservative. The volume fluxes we use are 


(B.l) 


= I {{hv}hu,m) 

d{hv}h,U,m) - b 

Similar to [5], a criterion for discrete entropy conservation is 

[gf (F™'+G™') = m + mf s = m+g {{mi m+g {{hv}} i&i, 

where the entropy potential (j) is defined as 

(f = q- (/+ 9)-{F+g) = ^gh'^u + ^gh'^v, 


(B.2) 


(B.3) 


where /, g are the physical fluxes (2.2), F, Q are the entropy fluxes (2.7) and we have used the consistent auxiliary 
source term discretisation 

( “ \ 
vstfli-l/ 


s \= 


(B.4) 


The jump in entropy variables is 


^glH + 51^1 - {{u}} M - {{?^}} H' 

Wi = I m 


(B.5) 
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We show (B.2| holds explicitly for the x—direction 


[g f i?-' = g {{hu}} [/rl - {{hu}} {{«}} M - {{hu}} {{.;}} H + {{hu}} {{«}} M 

+ g N - Ig {m} M + {{hu}} {{^}} H + g {{/r«}} [61 

= 9 {{^w}} 1^1 - ^9 {{^^}} M + 9 {{^}}^ N + 9 {{/im}} [61 

= 9 {{^w}} 1^1 - ^9 {{^^}} M + 9 {{^}} Ihuj - g {{h}} {{m}} [hi + g {{hu}} [61 

= 9 {{hu}} [hi - ^g {{h^}} [«! + ^[h^ul - g {{hu}} [hi - g {{h}} {{u}} [hi + g {{hu}} [61 ^ ^ 

= glhM - ^9 {{h^}} N - 9 {{/i}} {{^^1} [/il + 9 {{hu}} [61 
= glhM - IglhM + Ig {{«» Ih^l - 9 {{h}} {{«}} [hi + g {{hu}} [61 
= ^glh'^uj + g {{u}} [61 = l(l)j+g {{hu}} [61, 

and conclude that is entropy conserving. The y—direction is treated analogously to show that G’'°* is also 
entropy conserving. ■ 


Appendix B. 2. Proof of Property 3 


Proof. To verify the well-balancedness of the scheme we need to show that it solves the “lake at rest” problem 
correctly, meaning that the initial conditions h + b = const and u = u = 0 are preserved for all time. This happens 
if the discrete time derivatives vanish. It is immediately satisfied the discretised continuity equation since u or u 
are factors in all the terms. However, it is not immediately clear for the momentum equations. We show here in 
detail that the discrete time derivative vanishes for the hu equation, the hv equation is handled analogously. 


In the following proof we make extensive use of the Hadamard product notation (3.281 for the component-wise 
multiplication of matrices and define component-wise powers of nodal values by 


:= o W. 


(B.7) 


If we fully expand the flux-differencing form (4.2) of the hu equation by using (3.511 and the cubic forms from 
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we get the following scheme 


J(h o u)t + |(D(y^ o h o u^) + (-y^ oho u^)D^) 


+ -{y^o D(h o u^) + h o o D(y^) - h o o (yg)D^ - y^ o (h o u^)D^) 

+ -(h o u o D(yj; ou)+y^ohouo D(u) — h o u o (y^ o u)D’^ — y^ o h o u o (u)D'^) 

+ - (u o D(yj; ohou)+yj;Ouo D(h o u) — (yg oho — y^ o u o (h o u)D^) 

+ |(h o D(y^ o h) + y^ o h o D(h) - h o (yg o h)D'^ - y^ o h o (h)D^) 

+ -(—D(x^ o h o u o v) + (x^ o h o u o v)D'^) 

+ -(—Xjj o h o V o D(u) — h o V o D(x^ ou)+hovoxgo (u)D^ + h o v o (x^ o u)D^) 

+ -(—X,, o D(h ouov) — houovo D(x,j) + h o u o v o (x^)D^ + xg o (h o u o v)D^) 

+ - (—u o D(x^ ohov) — x^ouo D(h o v) + (xg o h o v o u)D'^ + x^ o u o (h o v)D'^) 

+ o (y^ o D(b) + D(y^ o b) - y^ o (b)D^ - (yg o b)D^) 

= -|m-i (y, o {{h}}^ o - I (-y^ o {{h}}^ o M,) + S (F* - F2) + (G^ - G2) S. 

(B.8) 

The source term discretisation is a key factor in obtaining a well-balanced scheme. In accordance with our 
approach for the flux terms, we have also discretised the source term as a quadratic split form. To account for 
possibly discontinuous bottom topographies on element interfaces we have introduced an additional interface part 
that vanishes for continuous topographies. The full source term discretisation is 


9 


|h o b;,; « |h o (y,, o D(b) -f D(y,, o b) - yj o (b)D^ - (yj o b)D^ 


+ |m- 


(yr, O {{h}}^ O IbJj) “ I (y« ° O N, 


M' 


(B.9) 


where we use notation (4.11 for the jump in bottom topography and the average water height. Using the definitions 
of the numerical and physical fluxes, the interface terms of the strong form DG discretisation are 

s(f;-F2) + (g;-G2)s = 

+ S (^y^ O ({{u}}^ o {{h}} -h I {{h^}}) -x^o {{u}} o {{v}} o {{h}} - y,, o (^h o V o -h -h x,, o h o u o v)) 

+ (-y« o ({{u}}^ o {{h}} -h I o {{u}} o {{v}} o {{h}} -h y^ o (h o -h - x^ o h o u o v) S. 

(B.IO) 

Under the “lake at rest” initial conditions, u = v = 0, many terms in the momentum equation (B .81 vanish. The 
remaining terms are 


J o (h o u)t -h ^{gh o D(y,, o h) -h y,, o gh o D(h) - gho (yj o h)D^ - y^ o o (h)D^) 
+ o (y,, o D(b) + D(y^ o b) - yj o (b)D^ - (yj o b)D^) 

= (y^ o {{h}}^ o IbJj) -I- S (y^ o ({{h^}} - )) 

- (-y? o {{h}}r, ° Wt)) -h (-yj ({{h^}} - h^)) S. 


(B.ll) 


We first consider the volume terms in ( |B. 11 ). For constant total water height h + b = const, the exactness of the 
derivative operator implies D(h-hb) = (h-Hb)D^ = 0. The remaining terms cancel due to the metric identities 
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(3.18) and we have 


|h o (D(y^ o (h + b)) + y,, o D(h + b) - (y^ o (h + b))D'^ - yj o (h + b)D^) 
= |h o ((h + b)D(y^) - (h + b) o (y4)D^) = 0. 


(B.12) 


Since we allow discontinuities in the bottom topography at element interfaces we must account for the jump in 
water height and we cannot guarantee that = h?. Instead we have at each interface 


m}-h-= 


hf + hi 


-hi = 


hi - hi 
2 


= ^ I^l> 


(B.13) 


where hi denotes the inner and ho the outer value. The sign of (B.131 depends upon which interface we consider. 


It is positive for the right and top interface and negative otherwise. In combination with the sign applied by the 
matrix, S = diag(^,..., — which gives a negative contribution for the top right component and a positive 
contribution for the bottom left component, we can restate these terms as 


S (y^ o = -M 1 (y,, o ({{h}}j o , 

(-y« ° ® ° ({{'^}}r, ° Mr,)) M-b 


With the reformulation (B.14I we show that the interface contributions are zero 


- (yn o {{h}}^ o Ib]^) - M-i (y,, o ({{h}}^ o IhJj)) 

- (-y« o {{h}}^ o |b],,) M-i - (-yj o ({{h}},^ o M"^ 

= - M-i (y,, o ({{h}}^ o |b + 

- (-y« o ({{h}}^ o [b + M-i = 0, 


(B.14) 


(B.15) 


since h + b = const, even on interfaces. 

We have shown that the discrete time derivative of the J{hu)t is zero. The treatment of the J{hv)t contribution 
is also zero by an analogous argument. Thus, the scheme is well-balanced. I 


Appendix C. Algorithmic Description of the ESDGSEM 


We outlined the reformulation of the curvilinear DGSEM in the flux differencing form in Sec. |3.4[ Now we 
provide specific details and restructure the algorithms of a standard DGSEM implementation to incorporate the 
entropy stable approximation. The flux differencing form for the curvilinear ESDGSEM (4.25) is self-contained, 
but it seems to be a daunting task to implement. We demonstrate, however, that with a few extra procedures and 
a slight restructuring of a standard DGSEM time derivative routine it is straightforward to implement the newly 
proposed entropy stable scheme. 


For this discussion we focus on the computation of the time derivative on a single spectral element. The global 
time derivative is assembled by looping over every element in a mesh nn. To make the discussion concrete we 
utilize the DGSEM structure outlined in Ghap. 8.4 of the book by Kopriva but the discussion readily extends 
to any standard DGSEM implementation. 


To begin, we introduce the notation used throughout this section. We store the computed solution on each 
element, scaled by the Jacobian, in the array {Jwij^n}^=ojl^‘^n=i where N is the polynomial order of the approx¬ 
imation and nEqn is the number of equations. The time derivative, also scaled by the Jacobian, is stored in the 
array {Each element stores array information about its mapping (the Jacobian, metric terms, 
etc.) in the geom object. We adopt the notation of a period to denote access to a component of an object. We store 
the bottom topography contributions in the volume and on the boundary of an element in separate places (for con¬ 
venience). The array stores the bottom topography evaluated at the Legendre-Gauss-Lobatto nodes. 

The array {dbi^j^n\f=Q’P=o^n=i stores the volume contributions of the source term and the array {lblijD}^=o id=i 
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stores the jump in the bottom topography along each edge of a quadrilateral element. To compute the source term 
at each edge we also store the average of the computed water height in the array id=i- 

First, we outline the details of the source term discretisation. We divide the computation into two parts: surface 
and volume contributions. For the surface contributions we alter the routine Alg. 137 (EdgeFluxes) from [17] . This 
is done out of convenience because the EdgeFluxes procedure already has access to local information about an edge, 
its local ID, and the elements that border an edge. Because we are on an unstructured mesh, care must be taken 
when computing the jump in the bottom topography term. We take the element to the left of an edge, ei, to be 
the “interior” and the element on the right, 62 , to be the “exterior” so that definition of the | 6 ] terms is clear. Also, 
we denote the local side ID of an edge on each element by Si and S 2 respectively. We fill a temporary array from 
the bottom topography term depends on the local side index of the edge on the elements ei and 

62 , e.g., if Si = 1 then = ei.{ 6 i^o}iLo S 2 = 3 then bn — e2-{bi^N}iLo- Then, after the normal numerical flux 
at an edge is computed, one adds to the existing algorithm: 

ei-{lbh,8Af=o = bR-bL, 

62 - = ei- > 

Alg. 0 {BottomContributions) is a straightforward implementation of ( |4.5[ ) to compute the volume contributions 
of the source term. We note that one could write the volumetric computation of the source term under BLAS3 
architecture standards. Because the bottom topography does not depend on time we precompute and store the 
quantities later use. 


Algorithm 1: (BottomContributions) Computation of the volumetric source term contributions. 
Procedure BottomContributions 

Input: 11 polynomial derivative matrix 

r JL N,N,nEqn , n 

{d 6 ij.n}iJo,j= 0 ^=l ^ 0 

for i = 0 to A do 
for j = 0 to A do 

{surnmjm^l <- 0 
for fc = 0 to A do 

sumi ■<— sumi + Di^k * bkj 
sum2 sum2 + 

sums ■«- sums + Di^k * {{geom.yr,)k,j * bkj) 
sumA sumA + {(geom.y^)i^k * 
sums sums + Di^k * {{geom.Xr,)k,j * bk,j) 
sums •«— sums + {(geom.x^)i,k * * Dkj 

dbij^2 {geom.yn)i,j * sums — (geom.y^)ij * sum2 + sums — sumA 
dbij^s < - {geom.Xn)i,i * sums + {geom.x^)ij * sums — sums + sums 

Output: 

End Procedure BottomContributions 


We next present Algs. ^{HighOrder-xFluxDifference) a,nd^{HighOrder-yFluxDifference) needed to reformulate 
the volume contributions of the ESDGSEM into the subcell flux differencing form presented in (3.471 and (3.48). 
Routines to compute the volume fluxes and G"®*, given in (IX 


are straightforward to implement, so we omit 
an explicit algorithm. To simplify the indexing in each of the high-order flux differencing algorithms we assume 
that the procedures is passed an appropriate slice from the solution storage and metric term arrays. 

Now that we have outlined the source term discretisation and the high-order flux extensions we are prepared 
to present the main algorithm for the efficient ESDGSEM implementation. We restructure the routine Alg. 144 
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Algorithm 2: (HighOrder-xFluxDifference) Computation of the flux difference AF. 
Procedure HighOrder-xFluxDifference 

Input: 11 slice of solution on Gauss-Lobatto grid 

Input: {{yn)i}f^Q , {iXn)i}^^Q // metric terms on Gauss-Lobatto grid 
Input: Q II SBP matrix 


for i = 0 to A do 

for m = 0 to A do 

{{Vv}} 2 + iyv)m) 

^ I ((a:^)i -b K)m) 




nEqn 


n = l 
'i nEqn 

) n = l 

N N,nEqn 

dFi^n !• 

J i=0,n=l 


End Procedure HighOrder-xFluxDifference 


{MappedDG2DTimeDerivatwe) from [17| that computes the local time derivative on a curved quadrilateral element. 
We outline the explicit steps to change a standard DGSEM approximation to implement the ESDGSEM for the 
shallow water equations: 


1. Begin with Alg. 144 {MappedDG2DTimeDerivative) from [T7] that computes the local time derivative on an 
element. 


2. Remove the standard approach that computes the DG derivative (denoted Alg. 92 {SystemDGDerivative) in 

[H])- 


3. 


Insert the equivalent flux differencing formulation for the volume terms outlined in Sec. |3.4| and detailed in 
Algs. I^and^ 


4. 

5. 


Use the entropy stable numerical fluxes (4.231 and (4.24) at element interfaces. 


Multiply the precomputed parts of the volume source term discretisation by the water height and add the 
source term contributions at each element in the volume. 


6 . Multiply the precomputed parts of the surface source term discretisation by the average water height and add 
the source term contributions at each element edge. 

7. Alg. 1^ {ESDG2DTimeDerwative) summarises the reformulation of a standard DG derivative to the compu¬ 
tationally efficient flux difference form. 


To reiterate, we assume that the approximations to the bottom derivatives are precomputed and stored and the 
surface source term contributions are computed in an augmented EdgeFluxes procedure. 
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Algorithm 3: (HighOrder-yFluxDifference) Computation of the flux difference GA^. 


Procedure HighOrder-yFluxDifference 

Input: II slice of solution on Gauss-Lobatto grid 

Input: // metric terms on Gauss-Lobatto grid 

Input: Q II SBP matrix 


for ji = 0 to A do 

for m = 0 to A do 

{{ye}} ^ \ ((ye). + (ye)™) 
{{®e}} ^ I ((®e)j + (a;e)™) 


Output 


{ - "s n£L,qn 

dG,,^ ^ 

J n = l 

{ - \ nEqn 

dG,,„ \ +2Q,,„ (- {{y,}} * + {{*e}} * ({ip,.4:fr, {w/™.4:fr)) 

J n = l 

{ - >. N + l,nEqn 

dG j n f 

J j=0,n=l 


End Procedure HighOrder-yFluxDifference 







(c) h, N = 3 


(d) Potential Vorticity, N = 3 



(e) h, N = 7 


(f) Potential Vorticity, N = 7 


Figure 7: At left is a three dimensional visualization of the water height, h, at T = 2 where wave breaking has 
occurred and bores are generated. On the right is the approximate potential vorticity (5.51 where the color ranges 
between -0.02 (black) to 0.02 (white). Numerical artifacts near element boundaries are visible due to the local 
derivative computation, without post-processing, of (5.51. 
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Figure 8: Grid convergence study for the ESDGSEM approximation of the parabolic partial dam break configuration 
with At = 1/1500 at T = 1.5. The overlay of quadrilaterals represents the mesh and the thick black line represents 
the unbroken portion of the parabolic dam. The color ranges between 3 (blue) and 10 (red). 
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(a) T = 0.0 


(b) T = 0.5 




(c)r = 1.0 (d)r=1.5 

Figure 9: ESDGSEM approximation for the parabolic partial dam break at four times with N = 5 and At = 1/1500. 
The overlay of quadrilaterals represents the mesh and the thick black line represents the unbroken portion of the 
parabolic dam. The color ranges between 3 (blue) and 10 (red). 
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(a) T = 0.0 


(b) T = 0.5 


(c)T=1.0 (d)T = 1.5 

Figure 10: Three dimensional visualization of the ESDGSEM approximation for the parabolic partial dam break 
at various times with N = 5 and At = 1/1500. Here the interaction of the flow with the discontinuous bottom 
topography is clear. The z—axis of the plot is from 0 to 10. 
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Algorithm 4: {ESDG2DTimeDerivative) Efficient implementation of two dimensional ESDGSEM in curvi¬ 
linear coordinates. 


Procedure ESDG2DTimeDerivative 

T • S 7Tr/. . \ N ,nEqn / / Vnr T .3T1 


Input 

Input 

Input 

Input 

Input 

Input 


{i=o'n^=i II solution scaled by Jacobian 

r fn* es 1 N,nEqn^A r gg N,nEqn,4 

<r ’ ml , 1 Lr ’ ml // entropy stable numerical fluxes 

I J mO,n=l,/D = l’ t ^'= 0 , 71 = 1 ,/D=l 

geom // element geometry 
Q // SBP matrix 

// Gauss-Lobatto quadrature weights 

{dKj,n}^=ofJo!n=lAlbh,ID}’i^2t,ID=lA{{h}}^,ID}^J),ID=l H bottom contributions 


// ^—direction 
for y = 0 to A do 

^ {geom.y^)o,j * xFlux ({Wo,i,n}"fr) “ {geom.Xn)o,j * yFlux ({lEo.j>}"f ?") 

{Eiv.n}"f -s- {geom.yr,)N,j * xFlux ({Wiv,j,n}^f i") - {geom.Xr,)N,j * yFlux ({lEjv,i,n}^f?") 

{ - "s N + ^,nEqn 

dFi,„ > ^ HighOrder-xFluxDif ference ({Wi,j,n}^o^l" , Q, {(5eom.i/^)ij}^4o , {( 5 eom.a;,,)i,j}^p) 

J 2=0,n=l 

for i = 0 to do 

{ — >1 nEqn 

{J 1 Eoj, 2 } ^ {JlEo,i, 2 } -f ^igeom.yr,)o,j * Mj,i * {{/i}}j ,4 
{JlEoj.s} ^ {JlEo.i.s} - ^(fleom.a;^)oo * Mi,4 * {{/i}}i,4 
{JlEv.i,2} ^ {JlEv,i,2} -f ^{geom.y^)N,j * Mi ,2 * {{^}}i,2 

_ {JWN,j,3} ^ {^lEiV.i.s} - ^igeom.Xr,)N,j * Mi.2 * {{^}}i,2 

// 77 —direction 
for i = 0 to A do 

{Go,^)HT ^ -(fleom.j/5)i.o * xFlux ({lEao.nj^f?") + {geom.x^)i,o * yP/w* ({lEi,o,„}::fr) 

i- {geom.y^)i^N * xFlux ({Wi,]v,n}^f?") -f {geom.x^)i,N * yFlux ({lEi,Jv,n}^f 1 ") 


{ - 'l +l,nEqn 

dGj,„ \ i- HighOrder-yFluxDifference ({lEi,i,„}^o®^d\ , Q, {{geom.yifji^j}^^^ , {{geom.x^)ij}f^^ 

J j = 0 ,n=l V J > J J 

for J = 0 to do 

^ ^ * {de,,A 

L 7 L J n = l 

{jm.Anl^fr ^ ({cr.Csllfr - {gnAZT) 

{JlEi. 0 . 2 } ^ {JlEi,o, 2 } - ^{geom.y()i,o * M^.i * {{^}}i,i 
{JlEi. 0 , 3 } ^ {Jlii.o.s} + 2^(fleom.a;5)i,o * Mi.i * {{^}}i,i 
{JlEi.iv. 2 } ^ {Jlii.v, 2 } - :^{geom.y^)i,N * Mi,3 * {{^}}i ,3 
_ {JlEi.iV. 3 } ^ {JWi,V,3} + :^{geom.Xi)i,N * Mi,3 * {{h}}i ,3 

for i = 0 to A do 
for j = 0 to A do 

L ^ - {{jw.,q,n}:^r+g * * {d 6 ,.,. 4 :fr) 

Output: {JWij,„}^^ofJo'!n^i 


End Procedure ESDG2DTimeDerivative 
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